Load packages
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3 ✓ purrr 0.3.4
## ✓ tibble 3.1.0 ✓ dplyr 1.0.5
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
Check working directory
getwd()
## [1] "/Users/jillashey/Desktop/PutnamLab/Repositories/Astrangia_repo/Astrangia_repo"
Load data
# Load data
daily <- read.csv("data/DailyMeasurements.csv", header = T)
daily <- daily[-13,] # remove tank 14 because it cracked
daily$Timepoint <- as.Date(daily$Timepoint, "%m/%d/%Y")
daily$Tank <- as.factor(daily$Tank)
daily <- subset(daily, Treatment=="Heat" | Treatment=="Ambient")
Discrete pH calculations from Tris calibrations.
path <-("data/pH_Tris") #set path to calibration files
file.names<-list.files(path = path, pattern = "csv$") #list all the file names in the folder to get only get the csv files
pH.cals <- data.frame(matrix(NA, nrow=length(file.names), ncol=3, dimnames=list(file.names,c("Date", "Intercept", "Slope")))) #generate a 3 column dataframe with specific column names
for(i in 1:length(file.names)) { # for every file in list start at the first and run this following function
Calib.Data <-read.table(file.path(path,file.names[i]), header=TRUE, sep=",", na.string="NA", as.is=TRUE) #reads in the data files
file.names[i]
model <-lm(mVTris ~ TTris, data=Calib.Data) #runs a linear regression of mV as a function of temperature
coe <- coef(model) #extracts the coeffecients
summary(model)$r.squared #extracts the r squared
plot(Calib.Data$mVTris, Calib.Data$TTris) #plots the regression data
pH.cals[i,2:3] <- coe #inserts coefficients in the dataframe
pH.cals[i,1] <- substr(file.names[i],1,8) #stores the file name in the Date column
}
colnames(pH.cals) <- c("Calib.Date", "Intercept", "Slope") #rename columns
pH.cals #view data
## Calib.Date Intercept Slope
## 20210122.csv 20210122 -86.73848 1.0218902
## 20210126.csv 20210126 -86.23930 0.9322418
## 20210203.csv 20210203 -88.14064 1.0017943
## 20210205.csv 20210205 -88.23311 0.8416322
## 20210217.csv 20210217 -87.53924 1.1237563
## 20210221.csv 20210221 -85.87848 0.8968835
## 20210223.csv 20210223 -86.36221 0.9370372
## 20210225.csv 20210225 -85.86760 0.9234509
## 20210227.csv 20210227 -87.60427 1.0150785
## 20210308.csv 20210308 -89.42691 0.9977500
## 20210317.csv 20210317 -86.21007 0.9799588
## 20210325.csv 20210325 -87.82550 1.0238389
## 20210409.csv 20210409 -88.31599 1.1153461
## 20210411.csv 20210411 -53.62103 -0.1385311
## 20210420.csv 20210420 -68.46091 0.2421390
## 20210524.csv 20210524 -96.95526 1.1877897
## 20210527.csv 20210527 -98.21886 1.2378231
## 20210531.csv 20210531 -88.25076 0.8102227
#constants for use in pH calculation
R <- 8.31447215 #gas constant in J mol-1 K-1
F <- 96485.339924 #Faraday constant in coulombs mol-1
#merge with Seawater chemistry file
SW.chem <- merge(daily, pH.cals, by="Calib.Date")
Calculate total pH.
mvTris <- SW.chem$Temperature*SW.chem$Slope+SW.chem$Intercept #calculate the mV of the tris standard using the temperature mv relationships in the measured standard curves
STris<-35 #salinity of the Tris
phTris<- (11911.08-18.2499*STris-0.039336*STris^2)*(1/(SW.chem$Temperature+273.15))-366.27059+ 0.53993607*STris+0.00016329*STris^2+(64.52243-0.084041*STris)*log(SW.chem$Temperature+273.15)-0.11149858*(SW.chem$Temperature+273.15) #calculate the pH of the tris (Dickson A. G., Sabine C. L. and Christian J. R., SOP 6a)
SW.chem$pH.Total<-phTris+(mvTris/1000-SW.chem$pH.MV/1000)/(R*(SW.chem$Temperature+273.15)*log(10)/F) #calculate the pH on the total scale (Dickson A. G., Sabine C. L. and Christian J. R., SOP 6a)
Plot daily measurements
# By Treatment
temp.trt = ggplot(SW.chem, aes(x=Treatment, y=Temperature.C)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Temperature°C") +
theme(axis.text.x = element_text(angle = 90))
temp.trt
## Warning: Removed 9 rows containing non-finite values (stat_boxplot).
sal.trt = ggplot(SW.chem, aes(x=Treatment, y=Salinity.psu)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Salinity.psu") +
theme(axis.text.x = element_text(angle = 90))
sal.trt
## Warning: Removed 21 rows containing non-finite values (stat_boxplot).
pH.trt = ggplot(SW.chem, aes(x=Treatment, y=pH.Total)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("pH Total Scale") +
theme(axis.text.x = element_text(angle = 90))
pH.trt
## Warning: Removed 22 rows containing non-finite values (stat_boxplot).
light.trt = ggplot(SW.chem, aes(x=Treatment, y=Light)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Light µmol photons m^2 s^-1") +
theme(axis.text.x = element_text(angle = 90))
light.trt
## Warning: Removed 625 rows containing non-finite values (stat_boxplot).
flow.trt = ggplot(SW.chem, aes(x=Treatment, y=Flow.mL.sec)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Flow mL/sec") +
theme(axis.text.x = element_text(angle = 90))
flow.trt
## Warning: Removed 1040 rows containing non-finite values (stat_boxplot).
plot.trt <- grid.arrange(temp.trt, sal.trt, pH.trt, light.trt, flow.trt, ncol=3, nrow = 2, clip="off")
## Warning: Removed 9 rows containing non-finite values (stat_boxplot).
## Warning: Removed 21 rows containing non-finite values (stat_boxplot).
## Warning: Removed 22 rows containing non-finite values (stat_boxplot).
## Warning: Removed 625 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1040 rows containing non-finite values (stat_boxplot).
ggsave("output/Daily_Measurements.By.Treatment.pdf", plot.trt, width = 21, height = 21, units = c("in"))
# By Tank
temp.tank = ggplot(SW.chem, aes(x=Tank, y=Temperature.C)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Temperature°C") +
theme(axis.text.x = element_text(angle = 90))
temp.tank
## Warning: Removed 9 rows containing non-finite values (stat_boxplot).
sal.tank = ggplot(SW.chem, aes(x=Tank, y=Salinity.psu)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Salinity.psu") +
theme(axis.text.x = element_text(angle = 90))
sal.tank
## Warning: Removed 21 rows containing non-finite values (stat_boxplot).
pH.tank = ggplot(SW.chem, aes(x=Tank, y=pH.Total)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("pH Total Scale") +
theme(axis.text.x = element_text(angle = 90))
pH.tank
## Warning: Removed 22 rows containing non-finite values (stat_boxplot).
light.tank = ggplot(SW.chem, aes(x=Tank, y=Light)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Light µmol photons m^2 s^-1") +
theme(axis.text.x = element_text(angle = 90))
light.tank
## Warning: Removed 625 rows containing non-finite values (stat_boxplot).
flow.tank = ggplot(SW.chem, aes(x=Tank, y=Flow.mL.sec)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Flow mL/sec") +
theme(axis.text.x = element_text(angle = 90))
flow.tank
## Warning: Removed 1040 rows containing non-finite values (stat_boxplot).
plot.tank <- grid.arrange(temp.tank, sal.tank, pH.tank, light.tank, flow.tank, ncol=3, nrow = 2, clip="off")
## Warning: Removed 9 rows containing non-finite values (stat_boxplot).
## Warning: Removed 21 rows containing non-finite values (stat_boxplot).
## Warning: Removed 22 rows containing non-finite values (stat_boxplot).
## Warning: Removed 625 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1040 rows containing non-finite values (stat_boxplot).
ggsave("output/Daily_Measurements.By.Tank.pdf", plot.tank, width = 21, height = 21, units = c("in"))
test for differences between tanks and treatments
# By Treatment
mod.temp <- aov(Temperature.C ~ Treatment, data = SW.chem)
mod.temp
## Call:
## aov(formula = Temperature.C ~ Treatment, data = SW.chem)
##
## Terms:
## Treatment Residuals
## Sum of Squares 3375.159 8500.317
## Deg. of Freedom 1 1109
##
## Residual standard error: 2.768546
## Estimated effects may be unbalanced
## 9 observations deleted due to missingness
summary(mod.temp)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 3375 3375 440.3 <2e-16 ***
## Residuals 1109 8500 8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 9 observations deleted due to missingness
hist(mod.temp$residuals)
mod.sal <- aov(Salinity.psu ~ Treatment, data = SW.chem)
mod.sal
## Call:
## aov(formula = Salinity.psu ~ Treatment, data = SW.chem)
##
## Terms:
## Treatment Residuals
## Sum of Squares 18.7864 1163.1030
## Deg. of Freedom 1 1097
##
## Residual standard error: 1.029688
## Estimated effects may be unbalanced
## 21 observations deleted due to missingness
summary(mod.sal)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 18.8 18.79 17.72 2.77e-05 ***
## Residuals 1097 1163.1 1.06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 21 observations deleted due to missingness
hist(mod.sal$residuals)
mod.pH <- aov(pH.Total ~ Treatment, data = SW.chem)
mod.pH
## Call:
## aov(formula = pH.Total ~ Treatment, data = SW.chem)
##
## Terms:
## Treatment Residuals
## Sum of Squares 2.57166 32.70446
## Deg. of Freedom 1 1096
##
## Residual standard error: 0.1727421
## Estimated effects may be unbalanced
## 22 observations deleted due to missingness
summary(mod.pH)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 2.57 2.5717 86.18 <2e-16 ***
## Residuals 1096 32.70 0.0298
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 22 observations deleted due to missingness
hist(mod.pH$residuals)
mod.light <- aov(Light ~ Treatment, data = SW.chem)
mod.light
## Call:
## aov(formula = Light ~ Treatment, data = SW.chem)
##
## Terms:
## Treatment Residuals
## Sum of Squares 90.49 124593.05
## Deg. of Freedom 1 493
##
## Residual standard error: 15.8973
## Estimated effects may be unbalanced
## 625 observations deleted due to missingness
summary(mod.light)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 90 90.49 0.358 0.55
## Residuals 493 124593 252.72
## 625 observations deleted due to missingness
hist(mod.light$residuals)
mod.flow <- aov(Flow.mL.sec ~ Treatment, data = SW.chem)
mod.flow
## Call:
## aov(formula = Flow.mL.sec ~ Treatment, data = SW.chem)
##
## Terms:
## Treatment Residuals
## Sum of Squares 3167.9 679148.5
## Deg. of Freedom 1 78
##
## Residual standard error: 93.31148
## Estimated effects may be unbalanced
## 1040 observations deleted due to missingness
summary(mod.flow)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 3168 3168 0.364 0.548
## Residuals 78 679148 8707
## 1040 observations deleted due to missingness
hist(mod.flow$residuals)
# By Tank
mod.temp <- aov(Temperature.C ~ Tank, data = SW.chem)
mod.temp
## Call:
## aov(formula = Temperature.C ~ Tank, data = SW.chem)
##
## Terms:
## Tank Residuals
## Sum of Squares 3384.931 8490.545
## Deg. of Freedom 15 1095
##
## Residual standard error: 2.784587
## Estimated effects may be unbalanced
## 9 observations deleted due to missingness
summary(mod.temp)
## Df Sum Sq Mean Sq F value Pr(>F)
## Tank 15 3385 225.66 29.1 <2e-16 ***
## Residuals 1095 8491 7.75
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 9 observations deleted due to missingness
hist(mod.temp$residuals)
mod.sal <- aov(Salinity.psu ~ Tank, data = SW.chem)
mod.sal
## Call:
## aov(formula = Salinity.psu ~ Tank, data = SW.chem)
##
## Terms:
## Tank Residuals
## Sum of Squares 23.8552 1158.0342
## Deg. of Freedom 15 1083
##
## Residual standard error: 1.034062
## Estimated effects may be unbalanced
## 21 observations deleted due to missingness
summary(mod.sal)
## Df Sum Sq Mean Sq F value Pr(>F)
## Tank 15 23.9 1.590 1.487 0.102
## Residuals 1083 1158.0 1.069
## 21 observations deleted due to missingness
hist(mod.sal$residuals)
mod.pH <- aov(pH.Total ~ Tank, data = SW.chem)
mod.pH
## Call:
## aov(formula = pH.Total ~ Tank, data = SW.chem)
##
## Terms:
## Tank Residuals
## Sum of Squares 2.71815 32.55798
## Deg. of Freedom 15 1082
##
## Residual standard error: 0.1734663
## Estimated effects may be unbalanced
## 22 observations deleted due to missingness
summary(mod.pH)
## Df Sum Sq Mean Sq F value Pr(>F)
## Tank 15 2.72 0.18121 6.022 3.19e-12 ***
## Residuals 1082 32.56 0.03009
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 22 observations deleted due to missingness
hist(mod.pH$residuals)
mod.light <- aov(Light ~ Tank, data = SW.chem)
mod.light
## Call:
## aov(formula = Light ~ Tank, data = SW.chem)
##
## Terms:
## Tank Residuals
## Sum of Squares 7884.24 116799.31
## Deg. of Freedom 15 479
##
## Residual standard error: 15.61537
## Estimated effects may be unbalanced
## 625 observations deleted due to missingness
summary(mod.light)
## Df Sum Sq Mean Sq F value Pr(>F)
## Tank 15 7884 525.6 2.156 0.00705 **
## Residuals 479 116799 243.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 625 observations deleted due to missingness
hist(mod.light$residuals)
mod.flow <- aov(Flow.mL.sec ~ Tank, data = SW.chem)
mod.flow
## Call:
## aov(formula = Flow.mL.sec ~ Tank, data = SW.chem)
##
## Terms:
## Tank Residuals
## Sum of Squares 51655.4 630661.0
## Deg. of Freedom 15 64
##
## Residual standard error: 99.26771
## Estimated effects may be unbalanced
## 1040 observations deleted due to missingness
summary(mod.flow)
## Df Sum Sq Mean Sq F value Pr(>F)
## Tank 15 51655 3444 0.349 0.987
## Residuals 64 630661 9854
## 1040 observations deleted due to missingness
hist(mod.flow$residuals)
Plot by date
# Plot and save
# By Treatment
timepoint.trt <- ggplot(SW.chem, aes(x=Timepoint, y=Temperature.C, group=Treatment)) +
geom_line(aes(color = Treatment), size = 1)
# By Tank
timepoint.tank <- ggplot(SW.chem, aes(x=Timepoint, y=Temperature.C)) +
geom_line(aes(color = Tank), size = 1) +
scale_fill_brewer("Dark2")
plot.timepoint <- grid.arrange(timepoint.tank, timepoint.trt, ncol=1, nrow = 2, clip="off")
ggsave("output/Daily_Measurements.By.Date.pdf", plot.tank, width = 21, height = 25, units = c("in"))